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SECTION  1 


INTRODUCTION 


The  development  of  light  weight,  high  strength  advanced  com- 
posite materials  affords  the  designer  numerous  possibilities  for  the 
production  of  improved  high-performance  aerospace  structures.  In 
order  to  fully  exploit  this  progress,  analysis  capabilities  must 
be  expanded  to  keep  pace  with  advances  in  both  materials  and 
fabrication  technology.  This  comment  is  particularly  true  in  the 
case  of  sandwich  composites.  Sandwich  construction  is  an  attractive 
alternative  to  other,  more  common  structural  configurations,  since 
substantial  benefit  in  terms  of  load-carrying  capacity  can  be 
obtained  under  little  or  no  weight  penalty.  However,  modes  of 
failure  which  are  unheard  of  in  simpler  types  of  construction  are 
possible  in  sandwich  structures  due  to  the  inherent  complexity  of 
the  geometry.  Design  and  analysis  methods  which  predict  these  modes 
of  failure  accurately  and  economically  are  therefore  a necessity. 

A particularly  important  consideration  in  the  effective 

utilization  of  sandwich  composite  construction  is  the  development 

of  nonlinear  analysis  capabilities.  It  is  a well-known  fact  that 

the  shear  flexibility  which  is  typical  of  sandwich  core 

materials  places  severe  restrictions  upon  the  range  of  validity 

. 1* 

of  linear  analysis  methods  . Linearization  has  been  shown  to 

be  particularly  restrictive  in  the  stability  analysis  of  sandwich 

2 

shells  having  significant  curvature  . 

It  is  clear  that  the  finite  element  method  is  the  most 
promising  approach  to  the  development  of  analysis  tools  which 
will  be  adequate  for  the  consideration  of  sandwich  composite 
structures  of  a very  practical  and  general  nature.  The  piecewise 
nature  of  finite  element  representations  permits  the  treatment 
of  complex  and  irregular  geometries  in  a straightforward,  unified 
manner.  The  consideration  of  very  general  loading  systems  and 
boundary  conditions  also  presents  no  particular  difficulty. 

* Numerical  superscripts  indicate  references  listed  at  the  end 
of  the  report. 
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The  present  report  outlines  a procedure  for  the  finite 
element  analysis  of  sandwich  structures  by  both  linear  and  non- 
linear methods.  A theory  of  sandwich  plates  is  presented,  and 
some  important  aspects  of  discretization  using  higher-order  (cubic) 
finite  element  representations  in  parametric  coordinate  space  are 
discussed.  Several  examples  illustrating  the  accuracy  and  flexi- 
bility of  the  numerical  analysis  are  given. 

Since  the  sandwich  structure  analysis  system  described  is 
not  yet  fully  developed,  many  limitations  remain  and  a number 
of  desirable  features  have  not  yet  been  implemented.  The  limi- 
tations of  the  existing  methodology  are  discussed  and  a number 
of  suggestions  for  further  development  are  mentioned.  It  is 
thought  that  the  analytical  approach  adopted  herein  represents  a 
useful  foundation  for  the  further  development  of  a reliable  and 
comprehensive  mathematical  model  for  the  analysis  of  practical 
structures  which  employ  sandwich  composite  construction. 

1 . 1 Description  of  the  Problem 

Sandwich  construction  is  a type  of  built-up  panel  configura- 
tion characterized  by  a number  of  thin,  high-modulus  layers  (face 
sheets  or  skins)  which  are  separated  and  stabilized  by  thicker,  low- 
modulus  layers  (cores) . Such  a panel  can  be  made  extremely  light- 
weight, while  offering  considerable  bending  rigidity  due  to  the 
separation  of  the  face  sheets.  The  function  of  the  sandwich  core  is 
analogous  to  that  of  the  shear  web  of  an  I-beam  section  in  coupling 
the  response  of  the  face  sheets  by  transferring  shear  stresses 
between  them. 

While  sandwich  composites  are  an  extremely  attractive 
concept  due  to  their  light  weight,  their  use  in  practical  struc- 
tures introduces  a number  of  possible  failure  modes  which  are  not 
encountered  in  other  structural  configurations.  Specifically, 
these  additional  failure  mechanisms  are  attributed  to: 

1.  weakness  of  the  core  layer (s)  in  comparison  with  the 
face  sheets,  both  in  shear  and  in  compression 
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2.  inadequacy  of  the  face  sheets  for  resisting  loads 
near  fasteners  and  inserts 

3.  lack  of  a continuous  interface  between  face  sheets  and 
core  in  cellular  (honeycomb)  core  constructions 

4.  the  possibility  of  defective  bonding  between  layers. 

Each  of  the  above  properties  give  rise  to  a particular  local  mode 
of  failure  which  is  unique  to  sandwich  type  construction.  Since 
each  of  these  types  of  failure  must  be  considered  both  in  formulating 
and  qualifying  a particular  design,  extensive  analysis  is  required. 
Furthermore,  a separate  checking  procedure  for  each  possible  type 
of  failure  may  not  always  be  sufficient,  particularly  when  the 
structural  component  in  question  is  to  undergo  large  deflections 
or  inelastic  behavior. 

Thus,  it  is  clear  that  sandwich  composite  materials  present 
a number  of  unique  problems  for  the  designer  or  analyst.  Proper 
treatment  of  these  problems  requires  a detailed  analysis  of  the 
sandwich  geometry  and  state  of  deformation,  to  which  the  finite 
element  approach  is  ideally  suited. 

1.2  Scope  of  the  Analysis 

The  present  analysis  addresses  the  problem  of  static  defor- 
mation of  flat  sandwich  panels.  Material  behavior  is  restricted 
to  be  linear  and  elastic,  although  geometric  nonlinearity  is 
considered.  The  analysis  is  thus  applicable  to  problems  of 
large  displacements,  elastic  buckling,  and  postbuckling . 

Sandwich  face  sheets  are  idealized  as  thin,  layered  com- 
posites obeying  the  Love-Kirchhof f assumptions.  Isotropic  and 
orthotropic  thin  plates  can  therefore  be  considered  as  special 
cases.  The  sandwich  is  of  the  anti-plane  type  (non-direct  stress 
carrying) , with  normal  deformations  considered.  The  finite  element 
discretization  is  performed  in  such  a way  that  multicore  panels, 

"half -sandwich"  constructions,  or  thin  laminated  composites  are 
easily  considered  by  stacking  elements  or  eliminating  individual 
layers  within  an  element  as  required. 
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External  loads  applied  to  the  sandwich  panel  may  consist  of 


concentrated  forces  and  distributed  pressure  or  body  forces.  Due  to 
the  numerical  integration  which  is  necessary  for  calculation  of  the 
element  properties,  the  specification  of  more  complex  distributed 
forces  (i.e.,  hydrostatic  or  sinusoidal  distributions)  is  also 
permitted. 

The  finite  element  discretization  presented  is  based  upon  a 
Hermite  bicubic  displacement  approximation  in  isoparametric 
coordinates.  Skewed  panels,  or  sandwich  having  arbitrary 
curvilinear  boundaries  can  therefore  be  considered.  The  method 
described  represents  the  first  development  of  an  isoparametric 
sandwich  finite  element  which  has  the  properties  of  both 
completeness  and  full  compatibility  of  displacements  as 
well  as  bending  slopes.  Considerations  in  the  parametric 
representation  of  finite  element  geometry  in  the  undeformed  state 
are  discussed  in  Section  3. 

\ j I 

The  research  reported  here  also  represents  one  of  the 
first  applications  of  the  parametric  bicubic  formulation 
to  nonlinear  problems.  A similar  finite  element  treatment  of  sandwich 
analysis  has  been  presented  by  Monforton^,  but  is  limited  to 

rectangular  shapes  in  Cartesian  or  cylindrical  coordinates.  Signifi- 

S r 

cant  gains  in  efficiency  have  also  been  made  in  the  present  work 
for  the  computation  of  nonlinear  strain  energy  and  energy  gradients 
for  use  in  numerical  solutions  by  direct  function  minimization  methods. 
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SECTION  2 

SANDWICH  PANEL  ANALYSIS 


A brief  description  of  the  sandwich  theory  which  is  used  as 
a basis  for  the  finite  element  analysis  is  given  in  this  section. 
The  geometry  of  deformation  for  a typical  sandwich  panel  is  dis- 
cussed, and  the  total  potential  energy  is  formulated  in  terms  of 

the  displacement  field  unknowns.  The  theory  presented  is  similar 

3 

to  that  outlined  by  Monforton  , but  is  extended  to  include 
normal  deformations  within  the  sandwich  core.  Large  displacements 
are  considered,  but  the  development  is  limited  to  small  elastic 
strains  and  moderate  rotations. 


2.1  Kinematics  of  Deformation 

Consider  a flat,  three-layer  sandwich  construction,  as 
shown  in  Figure  1.  The  two  face  sheets  are  taken  to  be  thin 
laminates,  each  composed  of  a number  of  orthotropic  layers  of 
arbitrary  orientation.  The  sandwich  core  is  a thicker,  relatively 
flexible  layer  having  constant  thickness  and  material  orientation. 


Each  face  sheet  is  considered  to  deform  according  to  the 
Love-Kirchhof f assumptions;  that  is,  linear  filaments  originally 
straight  and  normal  to  the  face  remain  so  after  deformation,  and 
undergo  no  extension.  In  terms  of  the  displacement  variables  at  a 
general  point  in  the  face  sheet,  this  assumption  implies  that 


Uf  - uf  - zfwf>x 
Vf  = vf  - zfwfy 

Wf  = wf  ; f = 1,2.  (2.1) 

Here  upper  case  letters  are  used  to  denote  the  displacement 
components  at  any  arbitrary  point,  and  lower  case  symbols  to 
represent  those  within  the  midplane  of  a face.  The  coordinate 
is  measured  upward  from  the  face  sheet  midplane. 

In  the  core  layer,  it  is  assumed  that  linear  filaments 
originally  straight  and  normal  to  the  layer  remain  straight,  but 


5 


need  not  be  normal  to  the  deformed  layer.  This  state  of 
deformation  is  pictured  in  Figure  2.  The  core  displacement  state 
therefore  takes  the  general  linear  form 

U = u + z <p 
c c c 

V = v + z ip 
ccc 


W = w + z Y 
c c CA 


(2.2) 


where  z^  is  measured  from  the  midplane  of  the  core. 

The  parameters  uc,  vc,  wc,  4> , \p , x can  be  eliminated  from 
Equation  2.2  by  enforcing  continuity  of  displacements  between 
the  faces  and  core.  Perfect  continuity  between  layers  is 
assumed,  so  that  the  possibility  of  debonding  failure  is  not 
considered.  At  the  lower  bond  line,  z1  = t.,/2  and  z = -t  /2 


so  that 


1 wl' 

t t 

U1  " 2~  Wl,x  = uc  " 2^ 


V1  ' 2“ 


t-,  t 

w.  = v 
l ,y  c 2 y 


Wn 


'1  c 

Similarly,  for  the  upper  bond  line. 


= W - y-X. 


(2.3) 


u. 


w. 


+ 2 w 

= u 

+ 

JS.  * 

2 z , x 

c 

2 

= V 

+ 

t 

2 z,y 

c 

2 

= wc 

+ 

t 

(2.4) 


Solving  Equations  2.3  and  2.4  for  the  core  displacement  parameters 
and  substituting  in  Equations  2.2,  one  obtains 


Uc  “ 2 (u2  + ui  + fc2W2,x  ~ fclWl,x) 

+ ^ lU2  - U1  + 2~v 2 ,x  + r-"!,*1 
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1 


R 


vc  = j(v2  + vi  + Y~  W2/X 


-iw 

2^1,  x 


+ _ (V2-V;L  + t2w2fX  + tlwl  x) 


Wc  = I (W2+W1)  + t5,  (W2‘W1) 

c 


(2.5) 


Thus,  the  displacement  components  of  any  point  within  the  sandwich 
faces  or  core  are  completely  determined  by  the  six  components  of 
displacements  within  the  faces:  u^,  v^,  w^;  f = 1,2. 

Calculation  of  strain,  stress  and  strain  energy  within  the 
core  layer  is  simplified  considerably  if  Equations  2.5  are  replaced 
by  approximate  conditions  of  continuity,  based  on  the  assumption 
that  the  face  sheets  are  very  thin  in  comparison  with  the  core. 

The  displacement  continuity  equations  are  then  applied  at  the 
midplanes  of  the  faces,  with  the  resulting  core  displacement 
expressions  given  by 

z 


Uc  = I (U2  + Ul>  + tT  (u2_ul) 

c 


Vc  = I (V2  + Vl>  + ir  (v2  - Vl} 

c 


Wc  2 (w2  + Wl}  + tc  (w2  - wx) 


(2.6) 


Equations  2.6  are  used  in  what  follows  to  enforce  disnlacement  con- 
tinuity of  the  sandwich  forces  and  core(s). 


The  strain-displacement  equations  of  a thin  face  sheet, 


valid  for  large  displacements  and  moderate  rotations,  are  given 


by  Novozhilov  : 

e 


xf  Uf , x + 2 wf,x  " zf,xx 


^ 1 2 

£ .p  = U-  + t Wc  - Z,W, 

yf  f,y  2 f,y  f f,.yy 


Yxyf  Uf,y  + Vf,x  + wf,xwf,y  _ 2zfwf,xy; 


(2.7) 


f = 1,2. 

For  the  sandwich  core,  the  strain  energy  arising  from  inplane 
deformations  is  assumed  to  be  negligible,  and  therefore  only  the 


J 
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transverse  shear  and  normal  strains  need  be  considered.  Further- 
more, it  is  assumed  that  the  strain  energy  of  the  core  is  adequately 
represented  by  a linear  relationship  between  strain  and  displace- 


ment 

. Thus , 

the  strains  within  the  core  layer  are  calculated 

from 

Y = 

'xzc 

k ‘W 

C 

+ I(W2,X+W1,X)  + 

^ (W2,x-"l,x» 

Y = 

' yzc 

Ec(V2'V1) 

+ I(W2,y+wl,Y>  + 

if  (u2,y"wl,y) 

e = 

zc 

1 (W2-Wl) 

(2.8) 

u 

c 

It  will  be  convenient  in  what  follows  to  refer  to  Equations 

2.7 

and  2 . 8 

in  matrix  form.  Letting 

{£}T  = 
f 

|>f  eyf 

YxyfJ  ! f - 1'2 

(2.9) 

and 

{e}T  = 
c 

jjxzf  Yyzf 

ezcJ  ' 

(2.10) 

the 

strains 

are  of  the 

form 

{ £ } f = 

{£°}f  + zf 

{<}f;  f = 1,2 

(2.11) 

and 

(£>c  - 

(e° }_  + z 
c c 

{<}  . 
c 

(2.12) 

2.2  Stress-Strain  Equations 

Each  of  the  sandwich  face  sheets  is  taken  to  consist  of  a 
number  of  thin  layers,  each  of  which  possesses  a stress-strain 
relationship  of  the  form 


(2.13) 


where  1,2  denote  the  principal  material  axes  of  the  layer.  Equation 
2.13  can  be  referred  to  the  structural  reference  axes  x,y  by  a 
linear  transformation  of  coordinates,  with  the  result 


Here  the  notation  refers  to  a stress  within  the  kth  individual 

layer  of  face  f.  By  substituting  for  the  strains  from  Equation 
2.11,  Equation  2.14  is  integrated  through  the  face  sheet  thickness 

5 

with  respect  to  the  weighting  factors  1 and  to  yield 


(A 


■V2  „<k) 


‘ [i/2  Qijf  dV  £ “ X'2 


(2.17) 


The  coefficients  A^j  and  in  Equation  2.15  are  the  overall 
membrane  and  bending  stiffness  of  the  face,  respectively, 
including  extensional-shear  coupling  properties.  The  B. . 
represent  the  effects  of  coupling  between  membrane  and  bending 
action  due  to  asymmetry  about  the  midplane  of  the  laminate. 

For  a sandwich  core  with  orthotropic  properties,  the  stress- 
strain  relation  is  of  the  form 


1 T 

1 XZ 

) T 

! 

_ 

) yz 

/ 

I Z 

c 

G 0 
xz 


0 

0 


G 

yz 
0 E 


xz 


yz 


-*c 


(2.18) 


Using  Equation  2.12  to  evaluate  the  vector  of  strains,  and 
integrating  over  the  thickness  with  respect  to  weighting  factors 


1 and  zc  yields  the  following 


( Q) 

| 

'A  » B' 

| 

|E"j 

R, 

= 

— * mmm- 

B . D 

l" 

l J 

Ic 

1 

c ’ 

(2.19) 


•'  —TV.— 


[A]  = 


t /2 
/ ° 


-V2 


G 0 
xz 

0 


[B]  = 


t /2 
/ C 
-V2 


Lo 

G 


yz 


0 

0 


dz 


xz 


yz 


0 

0 

E 


z dz 
c c 


(2.21) 


iDl  = 


t /2 
/ ° 


-V2 


xz 


0 

G. 


yz 


0 

0 

E 


2 , 
z dz 
c c 


-»  c 


2.3  Potential  Energy  Formulation 

Since  the  problem  under  consideration  is  conservative,  it  is 


possible  to  deduce  a potential  energy  of  deformation,  ir^,  such 


that  the  necessary  conditions  for  equilibrium  are  defined  by  the 
£ U 

condition 


6tt  = 0 . 
P 


(2.22) 


In  particular,  the  potential  energy  is  given  by 


tt  = U - W 
P 


(2.23) 


where  U represents  the  strain  energy  stored  in  the  body  due  to 
its  deformation,  and  W is  the  potential  of  the  forces  applied  to 
the  body,  both  expressed  in  terms  of  displacement  functions. 

In  terms  of  Lagrangian  stress  and  strain  functions,  the  total 
potential  energy  has  the  form 


" = ///  \ U>T  (o}dV  - //  { u } T { P } d A , 

^ V 

o S° 

a 


(2.24) 


where  S°  is  the  portion  of  the  surface  over  which  loads  are 
prescribed,  and  the  extent  of  integration  is  the  initial 
(undeformed)  volume.  The  vector  {P}  consists  of  generalized 
external  forces  conjugate  to  the  components  {u}  of  generalized 
displacement. 
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For  the  present  case,  the  stress  and  strain  functions 
in  the  sandwich  face  sheets  are  defined  in  Equations  2.7,  2.11, 
and  2.15;  in  the  core,  they  are  given  by  Equations  2.8,  2.12, 
and  2.19.  Making  use  of  the  direct  stress  and  moment  resultant 
forms  of  the  stress  parameters,  the  total  potential  energy  of 
the  deformed  sandwich  panel  is  as  follows: 


- ff  (u}T{P}  dA,  (2.25) 

S° 

a 

where  the  strains  are  defined  by  Equations  2.7  and  2.8. 

Expanding  the  energy  functional  in  terms  of  displacements 

gives 


7 T 

P 


¥'  f-P  <w2  * "vv2  * ^ 

A°  L c c c 

c 

+ 'W'W  (wl,x+w2,x)  * GyZc(v2-vl)(wl.yw2,y) 
+ 3Gxzctc(w2,x  + Wl,x  ~ wl(xW2,x* 

+ Kxe^  (w2/y  + wl,y  ' wl,yW2/y)]  dA 

+ |£S  //  {Allfu2f,x  + A22fVf,y  + A66f(uf,y+Vf,x: 


+ 2Ai 2fuf , xvf ,y  + 2(A16fuf,x  + A26fvf,y)(uf,y+v£,x) 


-2 


BllfUf ,xWf , xx  + B22fvf,y  f ,yy+2B66fWf ,xy ^Uf ,y+vf ,x* 


+ B. (v.  +u,  w,  ) + (u,  +v,  ) 

12f  f,xx  f,x  f,yy  16f  f,xx  f,y  f,x 

+2B. . ,u,  w,  +B_,ew,  (u.  +v,  )+2B,,_u,  w, 

16f  f,x  f,xy  26f  f,yy  f,y  f,x  26f  f,y  f,xy 
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+Dllfwf,xx  + °22fwf,yy 


+ 4D---W? 

66f  f , xy 


+ 2D12fwf,xxwf,yy 


+ 4D16fwf ,xxwf ,xy  + 4D26fwf ,yywf ,xy 


+2A66fwf,xwf,y(uf,y+Vf 

+ A12f(uf,xwf,y  + vf ,ywf ,x  > + A16fwf,x  (uf ,y+vf .x* 


+ Allfuf,xwf,x  + A22fvf,ywf.y  +2A66fwf ,xwf ,y luf ,y  f,x 
2 . ..  „2  „ , « ..2 


+ 2A16fuf ,xwf ,xwf ,y  + A26fWf,y(Uf,y +V£,x> 

+2A26fvf ,ywf »xwf ,y 

* [Bllfwf,xxwf,x  + S22fwf,yywf,y  + 4B66f  wf ,xwf ,ywf ,xy 
+ B12f (wf,xxwf,y  + wf,yywf,x)  + 2B16f 'Wf ,xywf , x+wf , xxwf ,xwf ,y 


+ 2B26f(wf,xywf,y  + w^ 


'f,yywf,xwf,y>] 

+ I (Allfwf,x  + A22fwf,y)  +(A66f  + ¥l2f] > wf , A ,y 

+ (A16fwf,xwf,y  + A26fwf,xwf,y)}  " 

2 

-£  (pxuf  + pyvf  + P2Wfl  ^ 

1*1  no 


(2.26) 


Thus,  the  potential  energy  of  the  deformed  sandwich  panel 
is  expressed  as  a function  only  of  the  displacements  within  the 
face  sheets  and  their  spatial  derivatives.  The  displacements 
(and  therefore  the  potential  energy)  of  the  sandwich  core  are 
obtained  by  what  may  be  thought  of  as  a linear  interpolation 
between  the  face  sheets. 

The  potential  energy  as  formulated  in  Equation  2.26  is 
sufficiently  general  to  describe  a flat  sandwich  panel,  of  an 
arbitrary  shape,  undergoing  large  deflections,  which  may  be  due 
either  to  the  intensity  of  loading  or  to  buckling  instability. 
Dissimilar  face  sheets  are  considered  as  well  as  any  asymmetry  of 
either  face  sheet  about  its  own  midsurface.  It  should  be  noted 
that  in  obtaining  Equation  2.26  it  is  assumed  that  the  core  prop- 
erties G , G , E „ are  uniform  through  the  thickness  of  the 
xzc  yzc  zc  ^ 

sandwich.  This  corresponds  to  setting  [B)c  equal  to  zero  in 
i Equation  2.19. 


SECTION  3 


FINITE  ELEMENT  DISCRETIZATION 

The  potential  energy  of  a flat  sandwich  panel  undergoing 
finite  displacements  has  been  formulated  in  Section  2.  In  the 
following,  the  details  of  a finite  element  discretization  of  the 
structure  are  considered.  Interpolation  of  the  element  displace- 
ment state  is  established  using  the  natural  (parametric)  coordinates 
of  a finite  element.  The  enforcement  of  continuity  of  both  dis- 
placements and  transverse  slopes  between  finite  elements  and  the 
representation  of  arbitrary  undeformed  geometries  are  also 
discussed . 

3.1  Interpolation  of  Element  Displacements 

The  choice  of  a method  of  interpolation  for  the  displace- 
ment variables  over  a single  element  is  of  fundamental  importance 
in  the  formulation  of  structural  finite  elements.  Not  only 
should  the  displacement  field  be  well-represented;  the  computation 
of  strain  and  stress  information  by  differentiation  of  the 
interpolation  formula  must  also  yield  acceptable  accuracy  with  a 
minimum  of  elements.  For  isoparametric  elements  where  the 
undeformed  geometry  is  also  represented  by  interpolation,  a 
careful  choice  of  basis  functions  is  essential  to  modeling 
accuracy  and  efficiency. 

Due  to  the  adoption  of  the  Love-Kirchhof f assumptions  for 
the  sandwich  face  sheets  in  the  present  analysis,  the  potential 
energy  functional  contains  second  derivatives  of  the  transverse 
displacements.  In  the  application  of  the  finite  element 

discretization,  it  is  therefore  required  that  the  approximate 

2* 

transverse  displacement  field  be  of  class  C on  the  interior 
of  a single  finite  element,  and  of  class  across  interelement 

7 

boundaries  . A suitable  displacement  approximation,  having 
extremely  good  convergence  properties,  has  been  introduced  by 

* A function  is  taken  to  be  of  class  Cn  if  it  possesses  at  least 
n-continuous , nonzero  derivatives  within  the  region  of  interest. 
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Bogner,  Fox  and  Schmit  . However,  that  formulation  is  restricted 
to  rectangular  boundaries.  In  the  present  analysis,  the  discretiza- 
tion is  performed  in  natural,  or  parametric,  coordinates,  so  that 
any  restrictions  on  the  undeformed  element  geometry  can  be 
eliminated.  It  will  be  shown  that  considerable  accuracy  can  be 
obtained  both  in  displacement  anc|  stress  prediction  by  this 
method  of  discretization. 

Consider  first  the  problem  of  interpolating  a function  F(£) 

in  one  dimension  such  that  continuity  of  both  the  function  and  its 

slope  are  everywhere  preserved.  By  requiring  the  interpolation  to 

reproduce  both  function  and  slope  exactly  at  each  sampling  point 

. 9 

one  obtains  the  first  order  Hermite  interpolation  formula  , 


F(S)  =^itH0j(?)Fj  + Hlj(^)F?j]' 


(3.1) 


where  the  ^ ( 0 are  cubic  polynomials  satisfying  the  conditions 


(3.2) 


where  are  the  interval  endpoints  (see  Figure  3) . 

For  the  purpose  of  obtaining  an  interpolation  formula  in  a 
non-rectangular  two-dimensional  region,  a new  set  of  coordinates 
can  be  defined  within  a unit  rectangular  region  which  is  then 
related  to  the  true  element  shape  by  suitable  transformation  of 
coordinates.  This  new  set  of  "parametric"  axes  may  be  thought  of 
as  a set  of  curvilinear  coordinates,  imbedded  within  the  original 
geometric  region  (Figure  4).  A function  F(£,n)»  defined  on  the 
rectangular  region  -l£.(£  ,n)^l,  can  be  represented  by  the  inter- 

Q 

polation  formula  , 

F(5'r’>  [Hoi(S,Hoj("IFij  + Hli(5'Hoj('>)F5ij 

+ Hoi<5'Hlj,'')Fnij  + (3.3, 
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obtained  from  the  product  of  two  one-dimensional  formulas  having 
the  form  of  Equation  3.1.  On  the  interval  (£^=-1 , £2=+l)  , the 
Hermite  polynomials  have  the  following  forms: 


Hoi(5)  = I(1+^i^)2  (2-q5> 


Hli(C)  = Ci(l+5i5)2(l-5i5) ; i=l,2. 


(3.4) 


It  is  important  to  note  that  the  interpolation  takes  place  within 
the  rectangular  region  in  £,n  coordinates,  and  that  the  nodal 


parameters  F^ , F^ 


represent  partial  derivatives  with 


respect  to  the  natural  coordinates  defined  within  a single  element, 


The  derivatives  of  each  of  the  displacement  components 
with  respect  to  the  natural  components  (£,n)  must  be  related  to 
derivatives  in  the  (x,y)  coordinates  in  order  to  properly  establish 
the  conditions  of  interelement  continuity.  By  the  chain  rule,  the 
transformation  is  of  the  form 


ru 

U'S 

PS 

s 1 

0 

0 

0 

ru,  ^ 

u'n 

<U'CC 

x,a 

Y'n  1 
y'Ul 

0 

0 

v-7  - 

0 

u, 

— y- 

u, 

XX 

u'nn 

x'nn 

y'nn, 

x'n 

2x,  y , 

ir  n 

U, 

yy 

ISnJ 

Lx'Cn 

y,£n  1 

SS 

y'£y'n 

x,_y,  +x, ry , 
V n K*  n 

u, 

V.  xyj 

(3.5) 


The  above  transformation  can  be  obtained  from  the  equations  of  the 
region  to  be  considered  (for  example,  the  polar  coordinate  trans- 


formation in  the  case  of  a circular  shape) , or  by  interpolations 
based  upon  nodal  coordinate  data^S  in  the  present  investigation. 


the  transformation  indicated  in  Equation  3.5  is  obtained  by 
approximating  the  spatial  coordinate  variables  in  terms  of  the 
natural  coordinates  of  an  element  in  the  same  form  as  Equation  3.3, 
and  simply  differentiating  the  interpolation  formula.  Some 
details  and  implications  of  this  process  are  indicated  in  Para- 
graph 3.2. 
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In  the  present  development,  all  three  displacement  components 
in  each  face  of  the  sandwich  panel  are  represented  by  the  bicubic 
expansion  indicated  in  Equation  3.3.  Although  this  order  of 
interpolation  is  greater  than  the  required  linear  formula  for  the 
insurface  displacement  components,  the  additional  degrees  of 
freedom  permit  the  option  of  enforcing  continuity  of  membrane 
strains  in  adjoining  elements.  The  resulting  stress  predictions 
are  correspondingly  more  accurate,  and  exhibit  only  minor  dis- 
continuities (due  to  local  bending)  between  elements. 

The  displacement  approximation  for  a typical  finite  element 
is  cast  into  a convenient  form  as  the  inner  product 


u-(C,n)  = {HU,nU  {U_ } ; f = 1,2 


where 


(H(5,n) } - |hq1  (C)hq1  (n) , 

Hll(C)H01(n) , 

H01(?)Hll(n) ' 

Hll(5)Hll(n) ' 

H01^)H02  (n)  ' 

Hn(5)H02(n), 

H01(5)H12(n) ' 

H11 (’) H12 (n) ' 

H02(?)H02(T1)' 

H12(5)H02(n) ' 

H02(5)H12(n) ' 

H12(5)H12(n) ' 

H02(5)H01(n)' 

H12(5)H01(n) ' 

H02(5)Hll(n) ' 

H12(5)Hll(n_|J 

(3.6) 


(3.7) 


{uf}  - I ufn,  uf5n'  uf nil ' uf£nil ' uf  12 ' uf 512 ' ufnl2 ' ufCnl2' 
1x16 

uf22'  uf 522 ' ufn22'  uf5n22'  uf 21 ' uf521'  ufn2l'  uf 5n2l|  • 


(3.8) 


The  remaining  displacement  functions  v^(£,n)  and  w^(^,n)  have  similar 
forms.  A displacement  vector  for  the  entire  sandwich  element  is 
then  assembled  in  the  form 


{ue>t  =Ku1}t,  {u2)t,  {v1}t,  tv2}T,  Cw1}T,  tw2>TJ  • 

1x96  L 


(3.9) 


The  resulting  sandwich  element  therefore  possesses  96  external  degrees 
of  freedom. 

It  should  be  noted  that  since  the  unknown  displacements  within 
the  entire  sandwich  panel  are  expressed  in  terms  of  the  face  sheet 
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displacements,  the  consideration  of  multicore  constructions  presents 
no  particular  problem.  Two  or  more  sandwich  elements  can  be  "stacked" 
to  represent  multicore  geometry,  simply  by  the  proper  specification 
of  element  connectivity  data.  The  displacement  vectors  (Equation  3.9) 
of  stacked  elements  are  joined  in  exactly  the  same  manner  as  for 
two  adjacent  planar  elements  during  the  accumulation  of  the  system 
equations . 

3.2  Parametric  Mapping  Considerations 

In  most  applications  of  the  isoparametric  element  formulation 
(e.g.,  solids,  plane-stress  elements),  only  the  physical  displace- 
ments are  needed  as  nodal  parameters.  However,  the  requirement  of 
slope  continuity  in  elements  based  upon  the  Kirchhoff  assumption 
necessitates  the  use  of  derivatives  of  the  displacements  as 
nodal  variables  to  fulfill  the  conditions  of  interelement  compati- 
bility. Since  the  displacement  derivatives  are  computed  in  parametric 
coordinates,  it  is  instructive  to  consider  the  constraints  placed 
upon  the  geometric  mapping  by  displacement  compatibility  conditions 
on  the  interelement  boundaries. 

For  the  enforcement  of  displacement  (but  not  slope)  continuity 
on  an  interelement  boundary,  examination  of  Equation  3.3  reveals 
that  the  displacement  at  any  point  on  the  boundary  depends  upon  the 
nodal  displacement  values  and  parametric  derivatives  tangent  to 
the  boundary,  evaluated  at  the  endpoints  of  the  edge.  Thus,  it  is 
sufficient  to  require  that  the  parametric  coordinate  tangent  to  the 
interelement  boundary  be  the  same  for  any  two  adjacent  elements  on 
their  common  edge.  Such  a requirement  presents  no  difficulty,  even 
for  mappings  of  a very  general  nature . 

The  need  for  slope  continuity  between  adjoining  finite  elements 
is  a somewhat  more  subtle  problem.  Consider  first  the  establishment 
of  slope  continuity  between  adjacent  rectangular  elements  (Figure  5) . 
It  is  easily  shown  using  Equation  3.2  that  a matching  of  the 
parametric  derivatives  and  cross-derivatives  at  common  vertices 
produces  continuity  of  transverse  slopes  in  parametric  space  along 
the  entire  boundary^.  Thus,  if  the  mapping  x=x  ( £,  n)  and  y=y(£,ri) 
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possesses  continuous  first  derivatives  on  the  common  boundaries, 

the  physical  slopes  (w,  and  w,  ) are  made  continuous  as  well 

x y 

whenever  w . w . . and  w . . are  matched  at  the  boundary  endpoints. 

C1]  mi  £ni] 

For  non-rectangular  elements  (Figure  6) , compatibility  of  the 
physical  slopes  is  obtained  under  similar  conditions.  However, 
since  the  transformation  from  parametric  to  physical  coordinates 
may  now  have  nonzero  second  derivatives,  one  must  examine  the 
conditions  under  which  equality  of  the  cross-derivatives  w~  . . 
at  the  common  nodes  is  admissible.  For  two  adjoining  elements  of 
identical  thickness.  Equation  3.5  dictates  that  a matching  of 
w^^j  at  common  vertices  is  permitted  only  if  the  transformation 
of  coordinates  possesses  continuous  first  derivatives  on  the 
interelement  boundaries  and  continuous  second  derivatives  at  the 
nodal  points. 

Under  the  conditions  stated  above,  equality  of  w^ j , w^j/ 
w^j  and  w£-nij  ensures  the  continuity  of  w,  w,^  and  w,^  (and  hence 
that  of  w,^  and  w,^)  along  the  common  edges  of  adjacent  elements. 

For  most  problems  having  reasonable  geometries,  no  particular 
problems  are  encountered  in  establishing  an  acceptable  mapping; 
in  the  most  general  case,  conditions  of  compatibility  can  be 
satisfied  by  the  enforcement  of  linear  constraints  between  nodal 
variables  involving  parametric  derivatives. 

In  the  present  development,  the  geometric  mapping  between 
the  physical  and  parametric  forms  of  a single  finite  element  is 
based  upon  the  same  interpolation  formula  as  the  displacement 
approximation : 


x(£,n)  = (HU,n)  } {X} 
yU,n)  = (H(c,n)  )T  {y} 


(3.10) 


The  vectors  {X}  and  {Y}  contain  values  of  x,y  and  their  derivatives 
with  respect  to  the  parametric  coordinates  £, n evaluated  at  the  nodes 
of  an  element.  The  geometric  mapping  given  by  Equation  3.10  is 
identical  to  the  well-known  "surface  patch"  representation  introduced 
by  Coons'^.  Palacol  and  Stanton^  have  adapted  the  bicubic  patch 
method  to  the  analysis  of  thin  orthotropic  plates  and  achieved 
excellent  results. 
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By  means  of  Equation  3.10,  the  problem  of  representing 
irregular  geometries  in  the  undeformed  state  is  reduced  to  the 
selection  of  an  acceptable  set  of  mapping  parameters  {X}  and  {Y} 
at  the  nodal  points.  For  many  shapes,  the  parameters  may  be  obtained 
analytically.  As  an  example  consider  the  circular  element  sector 
shown  in  Figure  7.  Taking  r(n)  and  0(C)  to  be  linear  functions. 


r = r + =-  (l+n)Ar 
o 2 


6 = | (1-?) A0 


(3.11) 

the  polar  coordinate  transformation  yields  x(5,n),  y(5,n)  explicitly 
as 


x = [rQ  + i(l+n)Ar]  cos^-d-5)  A0 
y = [rQ  + i(l+n)Ar]  sini(l-£)A0 
The  vector  {X},  for  example,  contains  the  entries 
<X}T  = l^i f XU'  xni'  x5nl'  x2 

For  the  annular  sector,  from  Equation  3.12, 


* f • • • f 


c5n4j 


(3.12) 


(3.13) 


{X}1  =|^r0  cosA0 
^■Ar  cosA0 


r A0sinA0 
2 o 

i-ArA0sinA0 

4 


(rQ+  r)cosA0,  i-(rQ+Ar)  A0sinA0 , 


jAr  cosA0 

r + Ar 
o 


r 


Ar 


Ar 


, 5rArA0sinA0 

4 

, 0 

, 0 
, 0 

. 0 


(3.14) 


As  a second  example,  consider  the  elliptical  element  shown  in 
Figure  8.  By  taking  the  element  edges  along  concentric  elliptical 
lines  and  lines  of  constant  angle  0,  the  following  description  of 
the  boundary  curves  of  an  element  is  obtained: 
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only  on  the  nodal  values  and  w^^)  on  the  common  edge 
presents  no  particular  problem,  it  is  required  that  the  slopes 


in  physical  coordinates,  w ,x  and  w,^,  be  everywhere  continuous 


on  the  boundary, 
found  by  setting 


The  required  conditions  of  continuity  are 


( w,  ] 

1 (w,  ) 

1 X 1 

X 

II 

) W,  I 

1 ) w,  ( 

l y ) 

' \ y) 

(3.18) 
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(3.19) 


For  the  simple  case  shown  in  Figure  9,  the  mapping  between 
the  global  and  natural  coordinates  can  be  chosen  as  bilinear, 

I 


x*  = a(l-£)/2 
I 


y = b(l-C)/2  + d(l+n)/2 


(3.20) 


and 


x11  = c(l+£)/2 
y11  = d (l+n)/2 


(3.21) 


Using  Equations  3.20  and  3.21  and  performing  the  arithmetic 
gives  the  constraint  equation 


(”'c)  J'o/ 

!".J  L o 


-c/a  -bc/adl  |w,j. 


II 


nc:i 


(3.22) 


which  must  be  satisfied  on  the  entire  common  boundary.  The 
second  equation,  which  reflects  the  fact  that  n is  identical  in 
both  elements  on  their  common  edge,  need  not  be  considered 
further.  The  remaining  equation, 


II  c,  I b I. 

w'^  a(w,£;  + dw,n) 


(3.23) 
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must  be  expressed  in  terms  of  the  nodal  parameters  at  the 
points  A,B  to  arrive  at  the  final  constraint  relations.  Using 
the  element  interpolation  functions  (Equation  3.6),  Equation  3.23 
can  be  expressed  as  a polynomial  in  n,  whose  coefficients  are  the 
nodal  variables  at  A and  B in  elements  I and  II.  By  requiring 
the  coefficient  of  each  power  of  n to  vanish,  one  arrives  at  the 
discrete  form  of  the  constraint  relations. 


0 -2 


3*  -3— 
d Jd 

-3*  3— 

d Jd 


2—  -2 
zd  ^ 


(3.24) 


Similar  constraints  to  ensure  interelement  compatibility  of 
displacement  and  transverse  slopes  can  be  evaluated  for  more 
general  element  shapes  and  orientations.  Equation  3.19  is 
again  used  as  a starting  point,  but  the  required  constraints 
can  become  considerably  more  complicated  than  Equation  3.24. 

In  such  cases  the  calculations  are  best  performed  numerically, 
in  terms  of  the  mapping  parameters  of  an  element. 

It  should  be  noted  that  in  most  cases,  an  acceptable 
parametric  mapping  can  be  constructed  even  in  the  presence  of 
highly  irregular  geometries.  For  example,  the  above  case  can 
be  considered  without  the  imposition  of  linear  constraints, 
even  though  the  physical  boundary  slopes  approaching  the  edge 
at  x = 0 are  double-valued.  Instead  of  the  bilinear  form, 
consider  the  mappings 


x(5,n)  = | (1-35  + 3£2  - £3) 

y(S,n)  = | [1-35  + 352  - 53  + 4^  (l-n)) 


(3.25) 


and 


x(£,n)  = | (1  + 3£  + 352  + C3) 

y(£,n>  - | (l  + n)  (3-26) 

for  elements  I and  II  respectively.  It  is  easily  verified  that 
the  above  parametric  mappings  are  twice  continuously  differ- 
entiable across  the  common  boundary,  which  is  more  than  adequate 
to  ensure  interelement  displacement  and  slope  compatibility. 


SECTION  4 

NUMERICAL  CONSIDERATIONS 

The  selection  of  the  best  and  most  efficient  numerical 
analysis  is  an  important  consideration  in  the  implementation  of 
any  finite  element  procedure.  In  particular,  the  calculation  of 
refined  element  representations  and  the  evaluation  of  nonlinear 
effects  must  be  done  efficiently  if  the  analysis  program  is  to 
be  cost  effective. 

In  this  section,  a number  of  aspects  of  the  process  of 
setting  up  and  solving  the  finite  element  equations  are  discussed. 
Many  of  the  details  involved  in  implementation  of  the  finite 
element  technique  are  widely  accepted1"^,  and  will  not  be  repeated 
here.  The  computational  techniques  considered  are  those  less 
common  in  finite  element  analysis,  or  developed  specifically  for 
use  with  the  present  method. 

4.1  Calculation  of  Element  Stiffness  Matrices 

The  computation  of  the  linear  stiffness  matrix  (corresponding 
to  quadratic  terms  in  the  potential  energy,  Equation  2.26)  is 
considered  below.  Evaluation  of  the  nonlinear  potential  energy 
and  the  energy  gradient  for  large-displacement  problems  is 
discussed  in  the  following  section. 

Since  the  element  geometry  is  not  predetermined,  it  is 
necessary  to  perform  numerical  integrations  to  obtain  the 
element  stiffness.  The  most  common  method  gives  for  the  stiff- 
ness matrix  the  integral14 

1 1 

[K]  = JJ[n]T[b]t[D] [B] [N] |j|d£dn  (4.1) 

-1-1 

where  [N]  is  a matrix  of  polynomial  functions,  [B]  is  the  strain- 
displacement  relationship  and  [D]  is  a matrix  of  elastic  constants. 
The  integration  can  be  carried  out  by  Gaussian  quadrature,  with 
the  integrand  being  evaluated  at  each  of  a number  of  sampling  points. 
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The  above  method  of  integration  is  not  practical  in  the  present 
case,  since  the  order  of  the  matrices  is  large  (the  element  matrix, 
[K] e , is  of  order  96),  and  because  many  integration  points  are 
required  for  the  exact  evaluation  of  the  higher-order  polynomial 
functions.  By  taking  advantage  of  repetitive  patterns  occurring 
in  the  polynomial  terms  in  the  original  functional  (Equation  2.26), 
the  number  of  operations  required  for  the  evaluation  of 
Equation  4.1  can  be  reduced  by  a factor  of  approximately  250, 
thereby  reducing  substantially  the  amount  of  computing  time 
required  for  setting  up  the  system  of  equations  to  be  solved. 

In  the  potential  energy  (Equation  2.26),  each  quadratic  term 
can  be  expressed  in  the  form 

T - if  ° u£,ij  U£,M  d*  dy'  (4’2> 

where  c represerfts--a  jnaterial  constant,  u^  and  u^  denote  the 
components  of  displacement  'tr^T-^^ , or  wf  within  a single  face 
sheet,  and  the  indices  i,  j.  k,  Z indicate  the  appropriate 
spatial  derivatives  of  the  displacement  components.  Restricting 
the  discussion  to  cases  for  which  c is  constant  over  the  element, 
the  use  of  Equation  3.6  gives 

T = | {u“}  TJ7<H,ijHH,]a}T  dx  dy  {u^} 

+ \ (u&  }T  fj  dx  dy  {u“}  (4.3) 

Thus,  regardless  of  the  values  of  a, 5,  only  the  16x16  matrix 

//  {H'ij }{H'ia}T  dx  dy  <4.4) 

and  its  transpose  must  be  evaluated  for  each  combination  of 
i,j,k  and  Z in  order  to  calculate  all  of  the  possible  combina- 
tions to  the  stiffness  matrix.  The  number  of  combinations  of 
i,j,k  and  Z is  further  reduced  when  only  the  upper  triangle  of 
the  element  matrix  is  considered,  since  the  component  matrices 
can  be  transposed  when  required  during  the  actual  assembly  of 
the  matrix. 
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The  numerical  integration  indicated  in  Equation  4.4  is  performed 
as  follows.  At  a particular  integration  point,  the  vector  {H(£,n)} 
(Equation  3.7)  and  its  derivatives  with  respect  to  the  coordinates 
(£,n)  are  computed,  and  transformed  by  the  chain  rule  to  (x-y)- 
derivatives.  The  resulting  vectors  are  multiplied  by  the  square 
root  of  the  product  of  the  Gaussian  weight  and  the  Jacobian 
determinant,  to  reduce  later  computations.  Finally,  all  possible 
outer  products  of  the  vectors  {H,. .}  and  {H,v .}  are  formed,  and  the 
sums  accumulated  over  an  nxn  grid  of  integration  points.  Having 
formed  the  component  matrices,  it  remains  only  to  multiply  through 
by  the  appropriate  elastic  constants  and  accumulate  the  products 
into  the  element  stiffness  matrix. 

It  is  significant  that  numerical  experiments  performed  using 
Equation  4.1  with  a 4x4  integration  grid  and  eliminating  the 
multiplication  of  all-zero  submatrices,  yield  computing  times  on 
the  CDC-6600  of  more  than  120  seconds  for  a single  element 
stiffness  matrix.  The  above  method  using  component  matrices 
requires  approximately  one-half  second  on  the  same  machine,  using 
the  default  level  of  compiler  optimization.  This  amount  of 
computation  is  approximately  the  same  as  that  expended  in  the 
stiffness  calculation  for  the  well-known  20-node  (60  d.o.f.) 
solid  element  using  a 3x3x3  integration  grid.  The  half-second 
computation  time  for  a single  element  matrix  can  be  made  even 
less  by  further  compiler  optimization  of  the  code. 

4.2  Evaluation  of  the  Nonlinear  Strain  Energy 

In  the  present  analysis,  the  solution  to  the  nonlinear  problem 
is  obtained  by  seeking  a minimum  of  the  potential  energy  functional 
directly.  Hence,  only  the  potential  energy  and  its  gradient  with 
respect  to  the  unknown  displacement  parameters  must  be  evaluated 


* Polynomial  terms  of  degree  2n-l  are  therefore  integrated  exactly. 
Integration  with  n=4  is  sufficient  for  most  geometries. 
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(see  Section  4.3).  Since  the  function  and  gradient  evaluation 
represents  a large  portion  of  the  total  solution  time,  and  must 
be  performed  a number  of  times,  it  is  important  to  organize  the 
computations  as  efficiently  as  possible.  The  numerical  inte- 
grations upon  which  the  calculation  of  the  linearized  stiffness 
matrix  is  based  (Section  4.1)  provide  a means  of  obtaining  all 
of  the  information  required  for  evaluation  of  the  nonlinear 
strain  energy  terms,  so  that  very  little  added  computation  is 
necessary. 

Consider  first  the  linearized  potential  energy,  which  can  be 
expressed  as 


J-{U}T[K]  {U}-{U}T{P} 

since  the  global  stiffness  matrix  [K]  is  symmetric, 
gradient  is  given  by 

Vir  (£)=  [K]  {U}-{P}  . 

P 


(4.5) 

the  required 


(4.6) 


The  direct  evaluation  of  Equations  4.5  and  4.6  is  the  most  efficient 
means  of  carrying  out  the  computation,  when  the  common  terms  and 
matrix  symmetry  are  taken  into  account.  Although  it  is  not 
necessary  to  form  a global  stiffness  matrix  (indeed,  this  feature 
has  often  been  cited  as  an  advantage  of  minimization  methods  of 
solution) , the  evaluation  of  the  potential  energy  and  its  gradient 
element-by-element  is  a relatively  inefficient  process.  A count 

it 

of  the  multiplications  and  input-output  operations  required  is 
evidence  of  this  fact. 


As  an  example  of  the  evaluation  of  the  nonlinear  terms  of 
the  strain  energy,  consider  the  cubic  term 


T C u?,ij  U£,M  uf,mn  dx  a* 


(4.7) 


* On  the  CDC  6000  series  computers,  where  memory-access  operations 
can  consume  somewhat  more  time  than  floating-point  arithmetic,  a 
similar  argument  holds  in  favor  of  the  present  method. 
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where  the  notation  is  that  used  in  the  previous  section.  For 

the  purpose  of  evaluating  the  term  itself,  it  is  convenient  to 
01  6 v 

evaluate  u_  . u^  , „ and  u’  directly,  using  the  vectors 
f,ij  f,KJo  f,mn 

{ H , i j } of  Equation  4.4  and  the  element  displacement  vector 
{U„ } (Equation  3.9).  This  computation  is  made  at  each  Gaussian 
quadrature  point,  and  the  weighted  sums  accumulated  to  evaluate 
T.  The  computational  expense  involved  is  small,  since  the 
vectors  {H,^}  have  already  been  calculated  at  each  integration 
point  during  the  linear  stiffness  matrix  evaluation.  Furthermore, 
only  a limited  number  of  terms  of  the  form  u^  . . appear  in  the 
energy,  and  repetitive  patterns  are  easily  taken  into  account. 


Evaluation  of  the  gradient  of  Equation  4.7  is  performed  in 
a similar  manner.  Since 


(4.8) 


the  cubic  term  can  be  rewritten  in  the  form 

T = cl[  up 

* C[t  u?,nm  dxdy(U«) 

+ C[l  u?,tj  dxdy{o|) 


(4.9) 


When  a is  different  from  $ and  Yr  the  gradient  of  T with  respect 
to  the  entries  of  {U^}  is  given  by 


3T  q 

= Cjl  (H,i:j)T  M uji]nn  dx  dy  (4.10) 

P Af  P 

and  again  the  integration  is  carried  out  numerically.  Since  each 
of  the  individual  terms  in  Equation  4.10  has  been  at  least 
partially  evaluated  in  evaluating  the  energy  function,  the  total 
gradient  evaluation  requires  only  a small  amount  of  additional 
computation . 
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Alternative  methods  of  organizing  the  calculation  of  nonlinear 

strain  energy  contributions  have  commonly  employed  very  large 

matrices,  whose  entries  are  associated  with  each  of  the  possible 

3 15 

combinations  of  nodal  variables  within  individual  elements  ' 

Such  techniques  often  require  excessive  amounts  of  storage  (for 
the  present  element,  a vector  of  length  9316  would  have  to  be 
saved  for  each  element  merely  for  the  evaluation  of  quartic 
terms  in  the  strain  energy)  and  input-output  time.  The  present 

method  uses  a minimal  amount  of  core  storage,  since  only  a single 

2 2 . 
vector  of  length  80n  is  nt^eded  for  each  element,  where  n is 

the  number  of  Gaussian  integration  points.  It  is  clear  that  the 

computational  procedure  described  here  represents  a very  effective 

approach  for  use  in  the  solution  of  problems  involving  geometric 

nonlinearities . 

4.3  Solution  of  the  System  of  Nonlinear  Equations 

In  the  present  analysis,  a solution  is  sought  by  direct 
minimization  of  the  discrete  potential  energy  with  respect  to 
the  undetermined  nodal  displacement  parameters.  The  solution 
process  is  therefore  one  of  unconstrained  function  minimization. 
Gradient  methods  of  minimization  have  been  shown  to  be  the  most 
powerful  class  of  solution  techniques  for  such  problems16.  Two 
such  methods  which  have  been  implemented  in  the  present  analysis 
are  described  in  the  following  section. 

4.3.1  Fletcher-Powell  Algorithm 

The  Fletcher-Powell  (or  variable  metric)  method  of 

minimization  bears  a close  resemblance  to  the  familiar  Newton- 

Raphson  iteration.  However,  the  Fletcher-Powell  technique  makes 

use  of  an  approximate  metric  in  place  of  the  matrix  of  exact 

second  derivatives  during  each  iteration.  The  algorithm, 

17 

originally  suggested  by  Davidon  and  improved  upon  by  Fletcher 
1 8 

and  Powell  , is  quadratically  convergent  and  possesses  extremely 
good  stability  properties. 
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Given  a function  F(X)  which  is  to  be  minimized,  the 
Fletcher-Powell  iteration  proceeds  as  follows: 


1.  An  initial  estimate  of  X_  of  the  solution  is 

o 


selected.  Usually  this  estimate  consists  of  a linear  solution, 
or  an  extrapolation  of  previous  nonlinear  solutions. 


2.  The  gradient  of  F is  computed,  and  an  initial 


search  direction  is  selected  along  the  direction  of  steepest 


descent;  that  is, 


S = - [H]  VF(X  ) 
o o o 


(4.11) 


where  [H]q  is  the  identity  matrix. 


3.  A value  of  a is  determined  in  such  a way  that 


F (X . + “S.)  is  minimized. 
1 1 » 


4 . The  vector  of  unknowns  is  updated  by 


-► 

X.^  = X.  + as. 

l+l  l l 


(4.12) 


5.  The  following  quantities  are  calculated: 


Yi  = vF(Xi+1)  - VF(Xi) 


(4.13) 


[M],  = a. 


■+•  >r 

s.  s. 

1 1 


i i -+  T 

S.  Y. 
l l 


(4.14) 


[N]i  =- 


([H]i  Y.)  ([H].  Y.) 


(4.15) 


Yi  [H]iYi 


[Hli+l  = Wi  + [M]i  + [Nli 


6.  A new  metric  is  computed  from 

(4.16) 

7.  A new  search  direction  is  determined  according  to 

Si+1  " -tHJi+1  VF(Xi+1)  (4.17) 

and  steps  3 through  7 are  repeated  until  convergence  is  achieved. 
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The  Fletcher-Powell  technique  exhibits  extremely  fast 
convergence  in  practice,  and  is  particularly  effective  when  a 
scaling  transformation  is  used  to  normalize  the  vector  X of 
unknowns'*'*’.  The  scaling  implemented  in  the  present  analysis 
adjusts  the  unknown  variables  according  to  the  diagonal  elements 
of  the  linearized  stiffness  matrix,  since  these  constitute  a rough 
approximation  to  the  second  derivatives  of  the  function. 

The  primary  disadvantage  of  the  Fletcher-Powell 
minimization  is  the  need  to  retain  a full  matrix  [H]^  at  each 
step,  since  storage  requirements  may  become  excessive  for 
larger  problems.  This  limitation  is  eliminated  altogether  with 
the  Fletcher-Reeves  algorithm,  as  discussed  in  the  following 
section.  For  smaller  problems,  however,  the  Fletcher-Powell 
scheme  is  possibly  the  most  powerful  of  the  gradient  methods 
of  minimization,  and  is  preferred  over  the  Fletcher-Reeves 
procedure . 

4.3.2  Fletcher-Reeves  Algorithm 

The  Fletcher-Reeves  (or  conjugate  gradient)  method 
. . .19 

of  minimization  is  similar  to  the  variable  metric  technique, 
but  utilizes  a less  involved  method  of  selecting  a direction 
for  search.  In  this  case,  the  information  retained  about  higher 
derivatives  of  the  function  consists  of  a single  vector,  rather 
than  a full  matrix.  The  current:  search  direction  S is  then 
selected  according  to 


|VP(Xi)  | 

|VF(^i_1) 


(4.18) 


= -VF (X . ) + 6.1.  . 

l li  l-l 


(4.19) 


replacing  steps  5 through  7 of  the  Fletcher-Powell  algorithm 
(Equations  4.13  through  4.17). 

It  can  be  seen  that  the  computer  storage  require- 
ments for  the  conjugate-gradient  algorithm  are  substantially 
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less  than  for  the  Fletcher-Powell  procedure.  However,  even 
though  the  Fletcher-Reeves  method  is  theoretically  quadratically 
convergent,  convergence  difficulties  can  arise  in  practice  due 
to  successive  roundoff  accumulated  in  the  £ vector  from  iteration 
to  iteration.  Periodic  restarts  of  the  iteration  are  therefore 
often  necessary,  and  convergence  may  require  substantially  more 
effort  than  is  theoretically  expected.  A thorough  discussion  of 
this  problem,  with  numerical  examples,  is  presented  in  Reference 
20. 

Substantial  improvement  in  the  performance  of  the 
Fletcher-Reeves  minimization  is  obtained  by  virtue  of  the  scaling 
transformation  mentioned  in  the  previous  section.  This 
modification  of  the  problem  serves  to  minimize  the  eccentricity 
of  the  function  space  under  consideration,  thus  improving  the 
numerical  behavior  of  the  iterations.  A number  of  examples 
demonstrating  the  effects  of  the  scaling  transformation  are 
cited  in  Reference  16. 


SECTION  5 


COMPUTER  PROGRAM 

The  sandwich  finite  element  analysis  reported  herein  has 
been  implemented  in  a computer  program,  which  is  briefly 
described  in  this  section.  Both  linear  and  nonlinear  analyses 
can  be  performed  with  the  program,  and  multiple-pass  capability 
is  provided  for  nonlinear  problems.  Although  the  code  has 
been  developed  on  the  CDC  6600  computer,  it  can  be  made  com- 
patible with  other  machines  with  relatively  few  modifications. 

5.1  Program  Size  and  Capacity 

The  sandwich  finite  element  code  presently  consists  of 
slightly  more  than  5000  FORTRAN  statements,  distributed  among 
68  subroutines  and  the  main  program.  Although  the  compiled 
program  occupies  less  than  70,000  (octal)  storage  locations 
on  the  CDC  6600,  linear  problems  having  up  to  6000  degrees  of 
freedom  and  nonlinear  problems  of  more  than  1000  degrees  of 
freedom  can  be  accommodated.  These  limits  are  easily  adjusted 
upward  by  modifying  appropriate  array-dimensioning  statements. 
Multiple  load  cases  may  be  considered  for  linear  analysis.  Three 
types  of  finite  elements  (plates,  cylindrical  shells,  flat 
sandwich  panels)  are  presently  available  in  the  program;  all  of 
these  are  based  upon  the  bicubic  interpolation  scheme  outlined 
in  Section  3.1. 

5.2  Program  Organization 

The  computer  program  is  arranged  in  modular  form,  and  utilizes 

21 

the  CDC  segmentation  loader  to  obtain  a maximum  of  flexibility 
in  ordering  the  computations.  The  segmented  structure  allows  for 
a greater  subdivision  of  the  computational  tasks  than  is  feasible 
with  the  conventional  overlay  organization,  and  hence  larger 
problems  may  be  considered  with  quite  modest  central  memory 
requirements.  Furthermore,  provision  for  creating  "copies" 
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of  certain  commonly-used  service  routines  makes  these  routines 
available  to  the  necessary  sections  of  the  program,  even  though 
they  are  not  retained  in  high-speed  memory  during  the  entire 
problem  solution. 

All  internal  data  is  normally  kept  on  disk  files,  in  contrast 
to  an  "overflow"  type  of  data  structure.  The  motivation  for  this 
is  the  large  number  of  degrees  of  freedom  involved  in  a single 
element  (96  for  the  flat  sandwich  panel  element) , so  that  problems 
considered  by  the  program  would  normally  be  solved  out-of-core 
regardless  of  the  data  organization.  The  handicap  for  small 
problems  is  minimal,  as  shown  in  Section  5.3. 

Five  different  equation-solving  routines  are  presently  avail- 
able to  the  finite  element  analysis.  For  linear  problems,  an 
in-core  band  solution  or  an  out-of-core  block  solver  are  provided; 
a conjugate  gradient  (Section  4.3.2)  algorithm  and  two  versions 
of  the  variable  metric  (Section  4.3.1)  solution  (in-core  and 
out-of-core)  are  implemented  for  nonlinear  analysis.  The 
method  of  equation  solution  is  chosen  according  to  problem  type, 
problem  size,  and  available  memory,  as  indicated  in  Table  1. 

The  present  computer  code  contains  partial  generation 
capabilities  for  nodal  coordinates,  element  connectivity  and 
boundary  condition  data,  to  reduce  the  effort  required  by  the 
user.  Geometry  data  and  distributed  loadings  may  alternatively 
be  described  in  equation  form  through  the  use  of  user-written 
subroutines.  A partial  restart  capability  is  also  provided  for 
use  in  nonlinear  analysis. 

5.3  Computing  Time  and  Memory  Requirements 

Memory  allocation  for  the  present  finite  element  analysis 
program  is  diagrammed  in  Figure  10,  according  to  blocks  of 
routines  which  represent  major  steps  in  the  computational 
procedure.  The  maximum  field  length  requirement  is  67740  octal 
words  of  memory,  with  this  being  determined  by  the  storage 
required  for  the  nonlinear  solution.  Linear  solutions  can  be 
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performed  using  only  62500  words  of  control  memory.  To  demonstrate 

the  efficiency  of  the  analysis  program,  the  code  has  been  tested 

on  a small  problem  and  compared  with  the  general  purpose  programs 

SAP-IV22  and  NASTRAN2 ^ . The  problem  is  that  of  a thin,  flat  plate 

loaded  at  the  center,  as  shown  in  Figure  11.  A four-element 

discretization  was  chosen,  resulting  in  eight  degrees  of  freedom 

for  the  SAP -IV  and  NASTRAN  analyses,  and  nine  freedoms  for  the 

present  method  (the  additional  degree  of  freedom  corresponds  to 

the  cross-derivative  nodal  parameter  at  the  interior  point) . All 

three  programs  were  executed  on  the  same  machine,  using  similar 

compiler  options.  Results  of  the  analyses,  including  computer 

run  times,  are  summarized  in  Table  2,  and  central  deflection 

results  are  compared  with  a closed  form  solution  given  by 
24 

Timoshenko.  Deflection  results  from  NASTRAN  and  the  present 
analysis  are  quite  good,  while  the  deflections  given  by  SAP-IV 
show  a significant  error.  It  can  be  seen  that  the  run  times 
required  for  the  program  here  compare  favorably  with  those  of  the 
other  two  programs.  It  should  be  pointed  out  that  the  sandwich 
computer  code  is  at  a slight  disadvantage  for  such  a small  problem, 
since  it  works  completely  out  of  core.  Furthermore,  stress 
computations  were  performed  at  25  points  in  each  finite  element  by 
the  sandwich  code,  while  stresses  were  obtained  only  at  the 
element  centroids  from  the  other  programs. 


* The  above  computer  memory  requirements  are  based  upon  capacities  of 
6000  degrees-of-freedom  for  the  linear  problem,  and  1000  for  nonlinear 
analysis,  on  a CDC  6600  machine. 


TABLE  2 

RESULTS  OF  COMPUTER  PROGRAM  COMPARISON: 
ANALYSIS  OF  THIN  PLATE  WITH  CONCENTRATED  LOAD 


Program 

NASTRAN 

SAP-IV 

Present 

Analysis 

Central  Memory 

140000g 

104000g 

64000g 

CP  Time 

10.1 

6.0 

6.2 

10  Time 

36.9 

14.7 

26.6 

PP  Time 

37.5 

21.9 

21.2 

Total  System  Time 

84.5 

42.6 

54.0 

Degrees  of  Freedom 

8 

8 

9 

Computed  Central 
Deflection 

0.2397 

0.2178 

0.2400 

Error 

2.0% 

11.0% 

1.9% 

i > - WWiJ  S**-* 


SECTION  6 

DEMONSTRATION  PROBLEMS 


Several  numerical  solutions  based  on  the  present  analysis  are 
presented  in  this  section.  Both  linear  and  nonlinear  problems 
are  considered.  Applications  of  the  method  to  multicore  sandwich 
constructions  and  nonrectangular  panels  are  illustrated.  It  is 
shown  that  good  accuracy  is  obtained  for  predictions  of  both  dis- 
placement and  stress  response,  even  in  the  presence  of  geometric 
singularities. 

6.1  Linear  Analyses 

A number  of  linear  analyses  are  presented  in  the  following 
sections.  The  bicubic  displacement  approximation  is  seen  to 
produce  very  smooth  stress  distributions  as  well  as  good  repre- 
sentation of  natural  boundary  conditions,  even  where  very  few 
elements  are  used  in  the  discretization. 

6.1.1  Multicore  Sandwich  Beam 

The  five-layer  cantilever  beam  shown  in  Figure  12  is 

considered.  All  three  faces  are  isotropic,  with  E = 10.0  x 10^ 

2 

lb. /in.  , u = 0.30,  and  are  0.025  inch  in  thickness.  Two  sets 
of  core  properties  (strong  and  weak  cores)  are  analyzed.  The  beam 
is  subjected  to  a transverse  end  load  of  10.  lb. 

The  finite  element  idealization  consists  of  eight 
sandwich  elements,  four  in  each  layer  of  the  beam.  The  use  of 
longitudinal  symmetry  and  the  enforcement  of  displacement  contin- 
uity in  the  central  layer  results  in  a total  of  237  independent 
degrees  of  freedom. 

Calculated  displacements  for  the  two  cases  considered 
are  listed  in  Table  3 . The  tip  displacements  for  the  weak-core 
case  indicate  the  importance  of  considering  normal  deformations 
in  the  core  layer;  end  deflections  of  the  upper  and  lower  faces 
are  different  by  more  than  four  percent  in  this  example.  Total 

51 


Figure  12.  Multicore  Sandwich  Bean 


TABLE  3 

FREE  END  DISPLACEMENTS  OF  MULTICORE  BEAMS 
A.  Weak  Core  Beam 


Upper  Face  (Unloaded) 

Central  Face 

Lower  Face  (Loaded) 


Displacement  on 

Displacement  on 

Centerline 

Free  Edge 

0.113 

0.106 

0.114 

0.106 

0.118 

0.108 

B.  Strong  Core  Beam 


Upper  Face  (Unloaded) 

Central  Face 

Lower  Face  (Loaded) 


Displacement  on 
Centerline 

Displacement  on 
Free  Edge 

0.0856 

0.0839 

0.0856 

0.0839 

0.0858 

0.0840 

stresses  in  the  lower  sandwich  face,  along  the  centerline  of  the 
beam,  are  shown  in  Figure  13.  Stress  distributions  for  both  beams 
are  quite  similar,  except  near  the  free  end  where  high  local  face- 
bending stresses  occur  for  the  weak-core  case.  At  the  clamped 
edge,  the  total  stresses  are  in  agreement  within  three  percent.  The 
total  moment  developed  at  the  root  section  is  computed  as  227.  in. -lb., 
based  on  average  membrane  stresses  in  the  two  outermost  faces . The 
actual  moment  due  to  the  tip  load  at  this  section  is  200.  in. -lb., 
so  that  the  natural  boundary  condition  is  well-represented  despite 
the  relative  coarseness  of  the  discretization.  It  should  be 
mentioned  that  the  above  stresses  have  been  computed  along  element 
edges,  and  that  the  values  given  for  x =0.,  5.,  10.,  and  20. 
are  true  nodal  stress  values,  which  are  notorious  in  most  finite 
element  formulations  for  their  sporadic  behavior.  The  primary 
reason  for  the  smoothness  of  the  stress  profiles  in  the  present 
analysis  is  the  enforcement  of  continuity  of  the  membrane  strains 
along  interelement  edges,  made  possible  by  the  use  of  derivative- 
degrees  of  freedom.  The  differences  between  adjacent  elements  in 
total  stresses,  evaluated  at  the  common  nodes,  are  0.15%,  0.80%, 
and  3.9%  for  the  strong-core  sandwich,  and  4.6%,  5.2%,  and  16% 
for  the  beam  with  a weaker  core. 

6.1.2  Skewed  Sandwich  Plate 

A flat,  rhombic  sandwich  plate  having  a skew  angle  of 

45  degrees  is  analyzed  for  response  to  a pressure  loading.  Such  a 

25 

problem  has  been  solved  previously  by  finite  elements  and  by  the 
2 6 

Ritz  method  . The  face  sheets  and  core  of  the  sandwich  have  the 
following  properties: 

Ef  = 1.  x 107  lb. /in. 2 vf  = 0.32  tf  = 0.025  in. 

G = 500.  lb. /in. 2 t = 1.00  in. 

c c 

The  transverse  compression  modulus  of  the  core  is  taken  to  be  very 

large,  so  that  normal  deformations  in  the  core  are  suppressed.  A 

2 

uniform  lateral  pressure  load  of  1.0  lb. /in.  is  applied  to  the 
entire  panel. 
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DISTANCE  ALONG  BEAM  AXIS  (in.) 


For  this  example,  the  entire  panel  is  modelled  by  sixteen 
sandwich  elements,  as  shown  in  Figure  14.  The  discretization  shown 
involved  328  independent  degrees  of  freedom. 

25 

Monforton  and  Michail  , using  a finite  element  approach 

developed  specifically  for  skewed  sandwich  panels,  obtain  for  the 

present  problem  a central  deflection  vjq  = 0.135  in.,  and  a maximum 

principal  moment  resultant  M = 25.3  in. -lb. /in.  Results  given  by 
2 6 

Kennedy  from  a Ritz  assumed-mode  approximation  are  wc  = 0.142 
and  M = 25.4.  Solution  of  the  problem  by  the  present  finite 
element  method  gives  w = 0.126  in.  for  the  maximum  transverse 
deflection,  and  M = 26.0  in. -lb. /in.  for  the  principal  moment 
(based  on  the  average  membrane  stress  in  each  face).  Thus,  the 
present  analysis  shows  the  panel  to  have  a slightly  greater  result- 
ant bending  stiffness  than  the  earlier  solutions. 

It  has  been  pointed  out  that  the  sandwich  finite  element 
developed  herein  allows  the  enforcement  of  continuity  of  membrane 
strains  between  adjacent  elements.  However,  the  local  face- 
bending stresses  depend  upon  the  curvatures  of  the  face  sheets, 
and  therefore  cannot  be  made  continuous.  It  is  interesting  to  note 
that,  for  all  four  elements  whose  corners  intersect  at  the  central 
node  of  the  skewed  panel,  corner  stresses  identical  to  four  decimal 
places  are  obtained  by  the  present  analysis,  even  though  the  local 
face-bending  stress  represents  approximately  25%  of  the  total 
stress  at  the  central  node  point. 

6.1.3  Axisymmetric  Annular  Plates 

As  an  example  of  the  accuracy  which  is  attainable  for 
nonrectangular  shapes  using  the  present  method,  a thin  (t=0.020  in.) 
annular  sheet  under  uniform  pressure  is  considered.  The  inner  and 
outer  radii  of  the  plate  are  10.0  in.  and  20.0  in.,  respectively 
(see  Figure  15) . On  the  inner  edge,  the  plate  is  fixed  to  a rigid 
shaft.  The  outer  edge  is  left  completely  free. 

The  geometrical  representation  for  this  problem  is 
described  in  Section  3.2,  and  the  natural  coordinates  ?,n  are 
shown  in  Figure  7.  Since  r = r(n)  and  0 = 0 (?)  in  this  case,  the 


are  zero  due  to  the  condition  of 
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nodal  parameters  and  w£nij  are  zero  due  to  the  condition  of 

zero  tangential  slope.  The  plate  is  modelled  with  a single  finite 
element,  resulting  in  a total  of  four  unconstrained  degrees  of 
freedom.  The  nonzero  nodal  parameters  correspond  to  w^j  and  w^j 
at  each  of  the  two  nodes  on  the  outer  radius. 


For  a total  load  of  1 lb.  (p  = 0.00106  psi) , the  one- 

element  solution  yields  a maximum  deflection  of  wq  = 0.1929 

inches,  and  a maximum  radial  stress  of  o . = 748.4  psi  at  the 

ri  24 

support.  An  analytical  solution  given  by  Timoshenko  shows 

that  w = 0.1989  and  o . = 1102.4  psi.  The  errors  in  transverse 
o n 

deflection  and  maximum  stress  are  therefore  3.0  percent  and  32.1 
percent,  respectively,  for  the  numerical  solution.  This  degree  of 
accuracy,  using  only  a single  element,  is  quite  good  particularly 
when  one  notes  that  the  average  aspect  ratio  of  the  finite  element 
used  is  slightly  greater  than  nine  to  one. 

The  same  problem  (with  r^  = 5.0  inches),  has  been 

solved  using  a 12  element  discretization  with  46  degrees  of  freedom. 

. 24 

For  this  model,  displacement  results  given  by  Timoshenko  are 

reproduced  exactly,  and  the  maximum  stress  at  the  clamped  support 

is  different  from  the  analytical  solution  b.y  1.4  percent.  The 

calculated  stresses  along  a radial  line  are  shown  in  Figure  16, 

normalized  with  respect  to  the  maximum  radial  stress.  As  in  the 

previous  problem,  a five  degree  sector  of  the  plate  is  considered 

in  the  model,  so  that  the  average  element  aspect  ratio  is  close 

to  unity. 

6.1.4  Circular  Sandwich  Panel 

A uniformly  loaded  circular  sandwich  plate,  shown  in 
Figure  17,  is  considered.  The  face  sheets  are  0.025  in.  aluminum, 

7 

with  = 10  and  v = 0.30.  Properties  of  the  core  are  G = 

£ 1 C 

26,000  psi  and  E = 10,  with  a thickness  of  0.450  inches, 
c z 

A previous  finite  element  solution  of  this  problem  has  been 
27 

obtained  by  Sharifi  , using  eight  axisymmetric  sandwich  elements 
in  the  radial  direction.  The  discretization  used  here  has  three 
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elements  on  the  radius,  with  the  innermost  of  these  containing  a 
geometric  singularity  where  the  element  edge  degenerates  to  a 
point . 


Stress  distributions  for  the  lower  face  of  the  circular 

sandwich,  shown  in  Figure  18,  agree  quite  well  with  the  results 

reported  in  Reference  27.  For  the  central  bending  moment,  the 

eight  element  discretization  of  Sharifi  yields  a value  of 

M - -340  in. -lb. /in.,  while  the  present  three-element  model  gives 

M = -360  in. -lb. /in.  At  radius  r = 20,  the  results  of  Reference 

27  are  Mr  - 485,  = 150;  the  present  sandwich  element  gives 

M = 417,  M.  = 155  (units  are  in. -lb. /in  in  each  case), 
r w 

This  particular  example  is  an  excellent  demonstration 
of  the  ability  of  the  Hermite  bicubic  element  to  maintain  good 
stress  accuracy  for  structures  having  curved  boundaries  as  well  as 
geometric  singularities.  It  should  be  noted  that  the  quality  of 
the  solution  does  not  diminish  in  the  vicinity  of  the  singularity, 
since  the  stress  computations  at  the  plate  center  are  just  as 
accurate  as  those  at  the  clamped  edge. 

6.2  Nonlinear  Analyses 

Examples  of  the  geometrically  nonlinear  analysis  of  sandwich 
plates  using  the  present  finite  element  method  are  described  in 
the  following  sections.  Load -deflection  data  are  presented  for 
the  case  of  a clamped  rectangular  panel,  for  which  a previous 
solution  is  available,  as  well  as  for  the  skewed  panel  whose 
linear  solution  has  been  given  in  Section  6.1.2. 

6.2.1  Rectangular  Sandwich  Plate 

A 50  inch  square  panel,  clamped  on  all  edges,  is 

subjected  to  a uniform  pressure  load  p.  The  sandwich  faces  are 

6 2 

aluminum  (Ef  = 10.5  x 10  lb. /in.  , uf  = 0.30)  and  0.015  in.  in 

thickness.  The  core  layer  is  1.0  in.  thick  and  has  shear  modulus 

G = 50000  lb. /in. 2 
c 

The  finite  element  model  for  this  problem  consists  of  four 


nonlinear  sandwich  elements,  as  shown  in  Figure  19.  Due  to  the 
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double  symmetry  of  the  problem,  only  one  quadrant  of  the  panel  is 
considered;  the  resulting  discretization  has  120  independent 
degrees  of  freedom.  Solutions  for  various  values  of  the  pressure 
loading  are  obtained  using  the  variable-metric  method  of  minimiza- 
tion (Section  4.3.1). 

The  load-deflection  path  for  the  present  problem  is 

shown  in  Figure  20  in  terms  of  a central  deflection  ratio  w /t 

3 2 2 c c 

and  the  loading  parameter  Q = 12  pa  (l-uf  ) /Ef tf tc  , where 
2a  = 50  in.  is  the  edge  dimension  of  the  panel.  These  results  are 
compared  with  previous  solutions  given  in  References  3 and  28.  The 
analysis  presented  in  Reference  28  involves  a perturbation  solution 
for  q in  terms  of  the  parameter  w /t  , and  is  limited  to  a two-term 

C C 

expansion  of  the  form  Q = C. (w/t)  + C_ (w  /t  ) where  the  ratio 

X c c zee 

w/t  is  on  the  order  of  one.  In  Reference  3,  a finite  element 
c c 

solution  is  given,  wherein  it  is  assumed  that  the  pressure  loading 
acts  at  the  midplane  of  the  sandwich;  such  an  assumption  is 
necessary  due  to  the  neglect  of  normal  deformations  in  the  sandwich 
core.  Antisymmetry  conditions  requiring  that  u^=  -u2 , v^=  -v^  can 
then  be  imposed  to  reduce  the  number  of  degrees  of  freedom.  It  is 
more  reasonable  to  suppose  that  the  loading  acts  over  one  face  of 
the  panel,  and  this  approach  is  taken  in  the  present  analysis. 
Figure  20  shows  that  the  present  method  yields  a solution  which 
is  slightly  more  flexible  than  that  of  the  previous  finite  element 
analysis,  and  substantially  more  so  than  the  perturbation  solution 
of  Reference  28. 

The  importance  of  a nonlinear  analysis  capability  is 
evident  from  the  load-deflection  path  for  the  present  example. 

The  deflection  predicted,  for  example,  at  Q = 30.  (p  = 27.7  psi) 
by  the  linear  analysis  is  more  than  55  percent  greater  than  that 
obtained  by  taking  geometric  nonlinearities  into  consideration. 

6.2.2  Skewed  Sandwich  Plate 

The  rhombic  sandwich  panel  considered  in  Section  6.1.2 
is  reanalyzed,  this  time  including  the  effects  of  geometric  non- 
linearities.  A sixteen  element  model  consisting  of  328  independent 
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degrees  of  freedom  is  employed  as  the  linear  analysis.  The 
Fletcher-Reeves  method  is  used  to  obtain  numerical  solutions. 


Load-deflection  results  for  the  skew  panel  are  shown 
in  Figure  21.  Numerical  solutions  at  pressures  of  1.0,  2.5  and 
5.0  psi  are  used  to  determine  the  curve,  with  the  starting 
estimate  for  each  solution  vector  extrapolated  using  all  previous 
solutions.  Although  the  Fletcher-Reeves  iteration  has  been  found 
to  perform  sluggishly  at  times  on  smaller  problems  (less  than  100 
degrees  of  freedom) , its  performance  in  the  present  problem  is 
quite  good.  Using  the  linear  solution  at  p = 1.0  psi  as  a start- 
ing point,  the  nonlinear  iteration  converges  in  only  17  iterations. 
The  solutions  for  2.5  psi  and  5.0  psi  are  converged  in  279  and  154 
iterations,  respectively,  using  a linear  and  then  quadratic 
extrapolation  for  the  starting  vectors. 

The  effect  of  nonlinearity  upon  the  response  of  the  skew 
panel  is  a great  deal  more  pronounced  than  for  the  square  plate 
considered  in  the  previous  example.  For  the  skew  panel,  the  central 
deflection  predicted  by  a linear  analysis  at  a pressure  of  only 
5.0  psi  is  more  than  67%  too  large. 

The  importance  of  the  computational  procedure  for  the 
element  potential  energy  and  its  gradient  (Section  4.2)  is  evident, 
since  more  than  two  such  evaluations,  on  the  average,  are  necessary 
at  each  iteration  in  the  minimization  solution.  For  the  present 
example,  the  average  execution  time  per  iteration  is  1.38  seconds; 
thus,  the  method  outlined  in  Section  4.2  is  seen  to  be  quite 
effective . 
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UNIFORM  PRESSURE  LOAD  (psi) 

Figure  21.  Central  Deflection  of  Skewed  Sandwich  Panel. 


SECTION  7 


SUMMARY  AND  RECOMMENDATIONS 

A parametric  finite  element  formulation  for  the  linear  and 
geometrically  nonlinear  structural  analysis  of  sandwich  composite 
panels  has  been  presented.  The  theoretical  development  includes 
effects  due  to  local  bending  deformation  in  the  sandwich  faces, 
and  both  transverse  shear  and  normal  strains  within  the  core 
layer.  Use  of  the  face  sheet  displacements  as  primary  unknowns 
in  the  discretization  permits  the  consideration  of  thin  laminates, 
multicore  sandwich,  and  "harf-sandwich"  constructions  as  well  as 
the  more  common  three-layer  geometry.  The  discrete  model  is 
based  upon  Hermite  bicubic  interpolation  in  parametric  coordinates, 
and  is  thus  applicable  to  panels  having  curvilinear  or  skewed 
boundaries.  Orthotropic  material  properties  and  effects  due  to 
membrane-bending  stiffness  coupling  are  also  accounted  for  in  the 
analysis.  Displacement  and  stress  results  obtained  using  the 
parametric  sandwich  finite  element  are  quite  good,  even  on 
relatively  coarse  meshes,  and  solution  accuracy  is  maintained 
even  when  modelling  curvilinear  shapes  with  degenerate  elements. 

The  results  summarized  in  this  report  represent  a preliminary 
state  of  development  of  a more  comprehensive  sandwich  composite 
analysis  capability.  The  discrete  formulation  described,  which  is 
capable  of  considering  a fairly  general  class  of  undeformed  panel 
geometries,  can  be  extended  to  model  sandwich  shell  geometry  as 
well  as  multiple  plate  or  shell  constructions  joined  at  arbitrary 
angles.  Provisions  for  including  stiffening  members  and  attach- 
ments to  surrounding  structure  are  important  as  well,  but  were  not 
implemented  in  the  present  effort.  Future  extensions  to  the 
present  analysis  should  also  include  consideration  of  material  non- 
linearity, and  dynamic  as  well  as  static  structural  response. 

The  research  described  herein  suggests  that  refined  finite 
element  approximations  such  as  the  parametric  bicubic  can  be 
successfully  adapted  to  the  general  nonlinear  analysis  of  structures 


of  a practical  nature.  Although  the  element  formulation  is 
relatively  complex,  the  resulting  numerical  performance  is 
extremely  good.  Computational  requirements  are  generally  modest 
due  to  the  high  solution  accuracy  which  is  attained  on  relatively 
coarse  grids. 
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